Dynamic subsurface changes on El Hierro and La Palma during volcanic unrest revealed by temporal variations in seismic anisotropy patterns

1 Introduction

Changes in volcanic systems are frequently accompanied by changes in subsurface stresses, often through the presence of melt pockets in the crust, inducing seismic anisotropy (e.g., Crampin, 1994; Gerst and Savage, 2004; Crampin, 1999; Illsley-Kemp et al., 2018). The measurement of shear-wave splitting provides one of the least ambiguous methods to probe subsurface anisotropy (Silver and Savage, 1994). This method takes advantage of the occurrence of two orthogonally polarised shear-waves, travelling at different speeds through an anisotropic medium, by measuring their time lag and polarisation (Silver and Chan, 1991). Individual results represent good measurements at a single location but are generally lacking vertical resolution (Savage, 1999). In studies using local events, varying results from earthquakes at different depths can reveal vertical changes of anisotropy, potentially caused by multiple anisotropic layers, which is often combined with a comparison to splitting from teleseismic events (e.g., Gledhill and Stuart, 1996; Pinero-Feliciangeli and Kendall, 2008; Schlaphorst et al., 2017; Schlaphorst et al., 2023).

Shear-wave splitting studies are mainly used to observe spatial variations in seismic anisotropy, but the use of this technique to look at temporal variations, such as those expected around volcanic eruptions, is not as common, as it requires an abundant wealth of data (e.g., Miller and Savage, 2001; Volti and Crampin, 2003; Gerst and Savage, 2004; Savage et al., 2015; Bianco et al., 2006; Roman et al., 2011; Baird et al., 2015). Whenever possible, the detailed study of temporal seismic anisotropy variations allows for identifying dynamic subsurface changes. In volcanic systems, the evaluation of shear-wave splitting measurement variations offers a helpful additional constraint that can be used in conjunction with other observations such as changes in seismicity and surface deformation, as well as geochemical analyses of erupted magmas and degassing events. Most volcanic eruptions are preceded by seismic precursors, which can have durations up to several years (e.g., Carracedo et al., 2022). This often provides a seismic source that can be used to track changes in seismic anisotropy and thus build a prolonged detailed temporal observation of changes in seismic anisotropy patterns, indicative of dynamic changes of the magmatic plumbing system (e.g., intrusive events), long before an eventual volcanic eruption takes place.

The westernmost islands of the Canaries, El Hierro and La Palma, are the sites of the most recent volcanic eruptions in the archipelago, posing a considerable risk to the local populations and, at the same time, representing a strong potential for high enthalpy geothermal resources (e.g., González, 2022; Serrano et al., 2023). Uplift of El Hierro in combination with local seismicity (González et al., 2013) is strong evidence of prolonged intrusions even after the volcanic eruption in 2011/2012 ceased (López et al., 2012; Benito-Saz et al., 2017). Even though not all instances of magmatic influx result in volcanic eruptions, the associated seismicity poses a seismic hazard and predictive models have generally high uncertainty. Therefore, a better understanding of continuing subsurface dynamic processes in the magma plumbing system of a volcano can help improve volcanic hazard analysis by decreasing the uncertainty of models, leading to improved mitigation of seismic and especially volcanic hazards. Moreover, a better understanding of the seismic response to magma transport and stalling within the critical zone may advance our forecasting capability of when shallow intrusions may evolve into eruptions.

In this study we present observations of temporal variations of seismic anisotropy beneath El Hierro and La Palma to constrain the structure and ongoing subsurface dynamics of the magmatic plumbing systems. We use local shear-wave splitting measurements from data recorded at five local seismic island stations in El Hierro and two in La Palma, and show that seismic anisotropy changed dynamically during the volcanic unrest that preceded, accompanied, and followed the recent eruptions on these two islands.

2 Setting

El Hierro and La Palma are the westernmost islands in the Canary Archipelago, a volcanic group of islands off the west coast of Africa, in the North Atlantic (Figure 1). The archipelago consists of seven major islands of volcanic origin, formed by the underlying Canary hotspot (e.g., Carracedo, 1999). P- and S-wave tomography shows a connection to a deeper reservoir in the mantle (Civiero et al., 2018; Civiero et al., 2019; Civiero et al., 2021), leading to volcanic activity. Due to the annual plate motion of the African plate (25 mm/yr in NE direction of 43°; Kreemer et al., 2014) over the hotspot, the emerging relative track results in a change of volcanic stages over the extent of the archipelago (Carracedo, 1999; Geldmacher et al., 2005; Gottsmann et al., 2008; González, 2023), which leads to a nearly-linear, age-progressive island chain towards the SW. After the formation of the next-oldest island to the east, La Gomera, the previously single-line volcanic belt appears to have been split into a dual-line oriented perpendicular to the general direction of the archipelago (Carracedo et al., 2001).

FIGURE 1. Map of the western Canary Islands of El Hierro and La Palma. (A) Overview of the archipelago. Bathymetry and topography are taken from ETOPO1 (Amante and Eakins, 2009), age contours (in Ma) are given by dark grey lines (Müller et al., 2008). (B) Location of the Canary Islands in the Atlantic Ocean. (C) Shaded Digital Elevation Model (DEM) of El Hierro, marking the location of the seismic broadband stations used in this study (white triangles) and their names, major rift (solid black lines) and the recent volcanic activity (2011–2012 eruption; red symbol). (D) La Palma shaded DEM model, marking the location of the Cumbre Vieja volcano, the Southern active volcanic rift in La Palma (2021 main eruptive fissure; red symbol).

The position of the islands relative to the hotspot also affects the crustal thickness, which increases from about 11.5 to 12.5 km in the west to about 20–30 km in the east, based on observations using Ps receiver functions (Martinez-Arevalo et al., 2013). This change is thought to be a product of cumulative intrusions and underplating as the islands get progressively older, as it happens in other volcanic archipelagos, such as the Cape Verdes (Lodge and Helffrich, 2006; Vinnick et al., 2012). The progressive age of the islands is also reflected in the topography of the islands throughout the archipelago, with the older eastern islands exhibiting razed morphologies and low elevations and the younger western islands exhibiting high, steep and youthful volcanic edifices. El Hierro and La Palma are the youngest islands with the most recent eruptions.

2.1 El Hierro

Volcanism on El Hierro has been active for as much as 1.2 Ma, first in a submarine environment followed by emergence at 1.1 Ma, with El Hierro becoming the youngest island of the archipelago (Carracedo et al., 2001). In July 2011, increased seismicity started in a central-western region onshore and moved offshore to the south in September 2011 and then to the north in October 2011. The last eruption took place from 10/10/2011 to 05/03/2012 offshore to the south of the island (Figure 1C) and marked the first volcanic eruption in the last 500 years. GPS measurements taken on this island show a period of uplift of about 27 cm until 2014, caused by a series of episodic magma intrusions during the post-eruptive period, ultimately being trapped at the Moho discontinuity at around 14–16 km (Klügel et al., 2005; González et al., 2013; Benito-Saz et al., 2017), and a stabilisation of the 2014 level afterwards (Figure 2D). However, seismicity is still ongoing and has shifted from swarms centred offshore to the north of the island to broader patterns onshore beneath the southern flanks, as well as offshore to the south and west of the island, hinting at ongoing subsurface activity. Two geodetically imaged reservoirs could be identified at 4.5 ± 2.0 km at the flank of the southern rift, and at 9.5 ± 4.0 km with a central source, of which sills are extending laterally, almost covering the whole island (González et al., 2013).

FIGURE 2. Seismicity maps during the observation periods in El Hierro [2010–2020; (A,B)] and La Palma [2021–2023; (C,D)]. The maps include all seismicity during the observation periods (A,C), as well as the selection of events that was used for the investigation based on magnitude and incidence angle (B,D). Seismic stations used in this study are indicated by triangles. Total number of events in each map is shown in the horizontal cross sections. Vertical movement from GPS measurements can be seen in the map insets for La Palma (B) and El Hierro (D). The apparent linear alignment of events at discrete depths is a result of rounded depth information in the catalogue, but does not affect our results.

2.2 La Palma

Volcanism on La Palma started at approximately 4 Ma in a submarine environment. After the island formation around 1.77 Ma the centre of eruption shifted from the northern part of the island to the southern Cumbre Vieja shield volcano, where all recent eruptions took place (Staudigel and Schmincke, 1984; Carracedo et al., 2001; Troll and Carracedo, 2016) and which marks the historically most volcanically active region of the Canary Islands (Carracedo et al., 2022). The island’s last eruption took place from 09/09/2021 to 13/12/2021 (Figure 1D), and ended up being the most voluminous in historical times at more than 200 million cubic metres (González, 2022; del Fresno et al., 2023). This eruption was preceded by seismic swarms as early as 2017, which marked the first time in over four decades that the area became active again (Torres-González et al., 2020). The seismicity is ongoing, though in contrast to El Hierro, the number of events with a magnitude over 2.5 is significantly lower. Here, GPS measurements show a slight deflation of about 4 cm during the eruption (Figure 2B; Charco et al., 2023). Three distinct magma storage levels could be observed around the southern part of the island at 9–13 km, a mushy intermediate one at 18–32 km, and at 33–38 km, based on geodetic and seismological results (del Fresno et al., 2023).

3 Methods

Subsurface seismic anisotropy can be investigated using observations from shear-wave splitting, as an initial shear-wave will split into two orthogonally polarised parts with different speeds when propagating through an anisotropic medium. The different speeds cause the two waves to separate in time and the polarisations and time lag will be retained even after leaving the anisotropic layer. Therefore, a measurement of the two splitting parameters, the fast polarisation direction (FPD, also abbreviated as ϕ) and the time lag (δ t), is used to characterise the anisotropy observed beneath a station. Similar to (Schlaphorst et al., 2017; Schlaphorst et al., 2023), we use the methodology by Silver and Chan (1991) to determine the splitting parameters, ensuring the stability of the solutions through cluster analysis of multiple analysis time windows (Teanby et al., 2004). In this method, a mathematical correction is applied to the data to remove the effects of splitting by a process of minimising the second eigenvalue of the covariance matrix.

In contrast to the more widely used measurements of teleseismic SKS phases that lose their initial source polarisation while travelling as P-waves through the outer core, the local S waves retain their source polarisation, which can complicate the detection of splitting parameters, as it is dependent on the focal mechanism of the earthquake. To increase confidence in the robustness of our results, we exclude all results with uncertainties larger than ±0.1 s (which approximates the 2σ uncertainty). Furthermore, surface reflections can distort measurements, the effect of which can be greatly reduced by limiting the selection of events to ones with steep incidence angles; the limit depends on the Poisson’s ratio and is approximately 35° in our setting (Evans, 1984; Savage, 1999). These influences in combination with often complex crustal structures often lead to small-scale changes of local anisotropic patterns (e.g., Huang et al., 2011).

In addition, the quality of the results is estimated following a set of criteria (Schlaphorst et al., 2017; Schlaphorst et al., 2023), which involves visual inspection. The criteria are: 1) the onset of the shear phase is clearly visible; 2) the signal-to-noise ratio (SNR) on the radial component is sufficient (SNR ≥ 3; SNR ≥ 10 for null results); 3) the corrected transverse component, shows a significant amplitude reduction, which is a sign of a successful reduction of the second eigenvalue of the covariance matrix that describes the particle motion; 4) the fast and slow shear waves show a similar shape; 5) the correction changes the particle motion from elliptical to linearised without cycle-skipping; 6) the error contour plot in the FPD and time lag plane shows a clear minimum of the second eigenvalue amplitude with a small well-constrained 95 percent confidence ellipsis; 7) both splitting parameters are stable over a suite of time windows with slightly changing start and end times; 8) a cluster around the best solution is clearly identifiable and few minor secondary clusters exist.

Due to significantly varying ray paths and different source polarisations of individual events, stacking results from local S phases to increase the SNR would bias the results. However, it is instead possible to sum results from different temporal or spatial bins in a circular histogram, also known as rose diagram. Results are weighted by their delay time and, in our case, grouped into azimuthal bins of 30°. The bin width is chosen to provide a good balance between azimuthal resolution and a useful amount of results per bin.

4 Data

On both El Hierro and La Palma, we use three-component seismic broadband stations from the IGN network (Figure 2; Supplementary Table S1) to ensure a large number of events with sufficiently steep ray paths. Due to noise levels, local events were chosen with a minimum magnitude of 2.5. However, the magnitude was lowered to 2.0 for La Palma data from 2022 to 2023 due to the low number of bigger events and lower number of successful splitting measurements. The data were filtered with a butterworth bandpass filter between 2 and 5 Hz to account for low SNR and oceanic noise. We tested for frequency dependency and found no significant variations in cases where signals could be observed with different boundaries, which has been observed in other oceanic settings (e.g., Di Leo et al. (2012); Schlaphorst et al., 2017).

4.1 El Hierro

In the case of El Hierro, an extensive catalogue exists, spanning 9 years of data since the eruption, from January 2011 to January 2020 (Figure 3A). We work with five seismic stations on the island. However, as most of the network was established after the onset of increased seismicity in July of 2011 and the start of the eruption on 10/10/2011, the data is insufficient to investigate events before and in the early stages of the eruption. Over the decade following the eruption, seismicity has shifted to the south and greater depths on average. Seismic swarms are clustered in a smaller area (lateral and vertical) until 2014 but have spread to larger extents in the following years. Due to the distance of the early seismic swarms to the then active stations there is only a small number of usable event-station pairs (Figure 2). Due to noise levels, local events were chosen with a minimum magnitude of 2.5.

FIGURE 3. Number of events and stations over time for La Palma (A) and El Hierro (B). Highlighted intervals mark subsets shown in Table 1 and subsequent figures. The subset naming convention consists of 1) island: EH (El Hierro), LP (La Palma); 2) year; 3) state of system and vertical GPS deformation: E (eruption), I (inflation), D (deflation), C (constant/no deformation). Horizontal stripes show the duration of the volcanic eruption. Usable events are station-event pairs within the 35° shear-wave window.

The data from El Hierro have already been presented in Schlaphorst et al. (2023). In total, out of 29,561 events, we find 2,263 events that fulfil the previously discussed selection criteria. Of those, we have data for 585 station-event pairs. Seismicity was first located around the central part of the island, before moving offshore to the south and the north during the eruption. In subsequent years, an increase in seismicity offshore to the west of the island could be observed. Here, we focus on the previously not discussed temporal variation of the dataset.

4.2 La Palma

In the case of La Palma, the network was established before the eruption on 19/09/2021. However, we do not yet have as many post-eruption years of data to add to the investigation (Figure 3B). In this study we use the two stations on the island that were made available to the scientific community, one of which (EHIG) is located directly above the main seismic swarm (Figure 2), allowing for a broad selection of event–station pairs.

The data from La Palma show a new set of results processed and assessed in a similar way to enable comparisons between the patterns observed on the two islands. In total, out of 11,461 events, we find 6,254 events that fulfil the selection criteria, of which we have data for 5,716 event-station pairs. In contrast to the El Hierro data, the general location of seismicity does not change during the observation time. We can identify two major source zones of similar lateral extent on top of each other in the southern half of the island and both provide a good event-station pairs. The main one (∼10,100 events) is located at a shallow depth between the surface and 20 km with hypocentre depths generally increasing towards the southeast of the region. The second largest (∼1,100 events) is deep, mostly located between 30 and 40 km and does not show any general depth variation depending on location. A third source zone of significantly sparser coverage (∼200 events) can be observed between the two major ones at an intermediate depth of 20–30 km. The median magnitude is below the minimum magnitude limit in the intermediate swarm (M_med = 1.9) but above in the shallow (M_med =2.5) and deep one (M_med =2.9).

5 Results

To emphasise temporal changes we have created time slices, which are grouped based on changes in seismic swarm location, GPS measurements, and/or splitting parameters chosen by visual inspection. A selection of individual time slices is shown in Figures 4, 5. Each slice shows a map with all individual measurements as bars representing the time lag by the length, and the FPD by the orientation of the bar. They are positioned on the midway between the event and the station and the paths are indicated by grey lines. In a further step, we summarise the results in circular histograms, called rose diagrams, which show general FPD. For El Hierro, the rose diagrams represent a combination of all five station results and for La Palma they are split between the two stations. A comparison of the ranges of event and splitting parameters between different time slices is shown in Figure 6 and the median values are listed in Table 1. The whole collection of shear wave splitting results are shown in Supplementary Table S2. In contrast to the often relatively low number of good quality results in teleseismic datasets (normally <5%; e.g., Schlaphorst et al., 2017; Schlaphorst et al., 2023), we find 215 results that pass the quality control (37%) in El Hierro and 908 (16%) in La Palma. The large number of results in the comparably short time span in La Palma can be explained by two factors: 1) the whole seismic network was active during the entire eruption when the largest number of events were recorded, whereas in El Hierro the network was not fully operational; 2) station EHIG is located directly above the main seismicity, accounting for 877 results. We do not find many null results, due to the fact that they are harder to identify in local seismic studies. This is due to 1) the retained source polarisation of the events and 2) the generally low SNR from small magnitude local events.

FIGURE 4. El Hierro splitting results for different time intervals (A–G). The results are plotted on the raypath at the midway point between the event and the station and summed up in rose diagrams, weighted by their time lags. Triangles depict stations active during the time span, red triangles show stations with good splitting results. Depth-time lag plots show the uncertainty in time lag. The colours of the events and splitting results are based on the depth of the event.

FIGURE 5. La Palma splitting results for different time intervals (A–H). The results are plotted on the raypath at the midway point between the event and the station and summed up per station in rose diagrams and weighted by their time lags. Triangles depict stations active during the time span, red triangles show stations with good splitting results. Depth-time lag plots show the uncertainty in time lag. The colours of the events and splitting results are based on the depth of the event. Rose diagram colours represent the shallow and deep source event swarms.

FIGURE 6. Parameter ranges of seismic swarms and groups for El Hierro (A–F) and La Palma (G–L). Details about median values can be found in Table 1 and about all shear wave splitting results in Supplementary Table S2.

TABLE 1. Details of event swarms and groups for both islands, including absolute number of days and events and median values of earthquake coordinates and splitting parameters. Note that for La Palma the groups are divided by station.

5.1 El Hierro

In El Hierro multiple significant changes in depth, FPD, and time lag, were observed over the entire decade since the eruption (Figures 4, 6, 7). General characteristics of individual swarms, such as their lateral location, event depth, as well as splitting parameters often tended to be similar within individual swarms (Figure 6; Table 1).

FIGURE 7. Temporal evolution of predominant FPD, vertical GPS variations and seismicity for El Hierro (A) and La Palma (B). Highlighted intervals mark subsets discussed throughout the paper and shown in Figures 3–6. Histograms show three colours for the number of all available events, usable events within the 35° shear-wave window, and good results used in this study. Rose diagrams for La Palma distinguish between events from the shallow (grey) and deep swarms (black).

A large number of events was observed between July 2011 and March 2012, before and during the volcanic eruption (∼12,700 events/757 events fulfilling all station-event pair selection criteria), accounting for about half the results recorded in total (97 results), as well as after the eruption during the time of the uplift recorded in the existing GPS network until January 2014 (∼9,000/1,400 events, 53 results). During that time, events with good results were confined to small distinct regions, the majority of which were located offshore north of the island, with a smaller extent of depth ranges; all except one were found around the base of the crust and/or at the upper mantle (15–26 km). Likewise, there is a strong uniform alignment of FPD within each individual swarm and, to a lesser degree, some distinct changes in time lag were found. A striking example of distinct variation in FPD alignment can be found between the swarms caused by two of the major magma intrusion events in 2012 and 2013 during the island uplift phase (Figures 4B, C; Figure 6E). Although the lateral distribution of sources is almost the same, the depth of the swarm shifted from 19–23 km to 15–19 km. At the same time, a ∼60° change of the FPD from a pattern aligned with a direction radially outward from the island centre to a pattern perpendicular to that was observed. This is accompanied by an overall decrease in time lag from a median value of 0.13 s, but reaching values as high as 0.43 s in EH12.I, to 0.11 s with a maximum of 0.25 s.

Between 2014 and 2018, seismicity was lower (∼3,000/53 events) and events exhibited generally lower magnitudes. Thus, only a small number of shear wave splitting measurements could be observed (6 results). During that time, no significant swarms were detectable, and the recorded results show no predominant FPD, while the general focus of epicentres shifted towards the south of the island with event depths of 13–30 km. The time lags generally increased with depth from 0.08 s at 13 km to 0.28 s at 24 km. This trend cannot decisively be confirmed due to the low number of results; however, a similar same trend with more results was observed during 2018, adding confidence to the observed pattern.

In 2018, the seismicity and number of results started to increase again (∼1,600/18 events, 13 results), without showing distinct swarms of seismicity. Events were located south of the island, mostly offshore, and were located at generally greater depths of 16–39 km. Time lags are similar with 0.08 s at 19 km to 0.30 s at 39 km, although larger time lags of up to 0.57 s can be found. In addition, a predominant FPD with N-S orientation was observed, though the direction is less well-defined than those in the seismic swarms prior to 2014.

Until August of 2019, a mostly similar pattern to the previous year continued. Events were located almost exclusively offshore in a cluster of smaller extent at the same depths. Although the numbers of events and results are slightly increased (∼1,900/14 events, 15 results) no major swarms could be observed. The time lags at shallower levels are increased (minimum of 0.13 s at 19 km) and slightly decreased at greater depths (maximum of 0.29 s at 36 km), while the FPD showed a weaker predominant pattern of NW-SE orientation.

The first major swarm of results in the time since the island uplift stopped at the end of 2014 can be observed from 16/09/2019 to 26/09/2019 (13 results). The events were located to the west and southwest of the island, again mostly offshore at depths of 30–36 km. The time lags showed a range independent of event depth but reaching high values of up to 0.41 s. The FPD shows strong alignment in NW-SE orientation.

A second swarm on 26/01/2020 (9 results) continued the general increase in seismicity. The epicentres were located mostly offshore south of the island at depths of 30–32 km. Time lags show moderate values of 0.15–0.23 s with one higher value of 0.58 s, while the FPD is uniformly aligned in N-S direction. The higher time lag value was recorded at the northernmost station, making it the only ray path with a majority located onshore crossing the island.

5.2 La Palma

In La Palma the majority of seismicity with a magnitude above 2.5 was observed from mid-September to mid-December of 2021 during the volcanic eruption and island subsidence (Figures 5–7; Table 1). Only few significant changes in event location and depth occurred during that time, the most notable being a slight increase in depth of the shallow swarm events between events before the eruption (LP21a.C; Figure 5A) and a few days after the eruption onset (LP21b.ED; Figure 5B), while the lateral location stayed the same within uncertainties. In addition, splitting parameters are mostly consistent during that time. After the eruption, the number of results is very small, so comparisons are not straightforward.

For station EHIG events are located at 7.5–16.0 km in the shallow swarm (842 results) and 31.2–37.5 km in the deep swarm (35 results). Due to its geographical location, results at station TBT only come from events in the deep swarm at depths of 28.0–38.9 km (31 results). Results from the shallow swarm show time lags that are usually <0.25 s (790 results), but there is one observation with a larger time lag of >0.5 s. At the deeper swarm we observed a slightly smaller portion of results with time lags <0.25 s (EHIG: 30 results; TBT: 25 results). The FPD for results at station EHIG is predominantly E-W, except for a 4-day time period in November, where a ∼40° change to an NW-SE orientation was observed (Figure 5C; Table 1). At station TBT, multiple patterns with strongly uniform orientation were observed (Figures 5A, C).

During 2022 and 2023, after the volcanic eruption, the number of results drops significantly, due to the fact that most events have magnitudes smaller than the 2.5 threshold. Even when lowering the minimum magnitude limit to 2.0, the number of results only increases from 4 to 25. The FPD at station EHIG shifted from predominantly E-W to NW-SE during a short interval of increased seismicity (LP22.C; Figure 5H) to N-S in 2023, although numbers are too small to draw definite conclusions of general orientation patterns.

6 Discussion

On both islands, results from temporal variations of shear wave splitting results are able to help explain the different stages of subsurface magmatic movement and the associated volcanic eruption. Especially on El Hierro where variations of event locations and splitting parameters in different time slices are significant, a complete summation of all results (see Schlaphorst et al., 2023) can only show parts of the detail and omits temporal variations entirely. Likewise, although less pronounced, some variations are visible on La Palma.

6.1 El Hierro

Evidence suggests that prior to the El Hierro eruption in May 2011, the deeper magma reservoir at the bottom of the crust and the uppermost mantle experienced a magmatic intrusion, leading to inflation, which resulted in the general uplift of the island (González et al., 2013). Based on elastic models of ground deformation derived from a combination of multifrequency and multisatellite data it is inferred that the magma chamber is located at a depth of around 9.5 ± 4.0 km, directly beneath the island and thus north of the location of the eruption. Changes in magmatic volume and resulting vertical crustal movement are likely to cause stress, which influences seismic anisotropy (e.g., Gerst and Savage, 2004). However, the first major increase of seismicity around the central island was unable to provide good event-station pairs due to the shallow incidence angle from large events (M > 2.5) to the active stations that provided good data at that time, and as such it was not possible to compute shear-wave anisotropy solutions.

Subsequently, starting in late September 2011, magmatic influx to a secondary smaller and shallower reservoir at around 4.5 ± 2.0 km offshore but close to the southern coast of the island led to an increase in seismicity to the south (González et al., 2013). Again, shallow incidence angles to active stations at the time prevented an analysis of robust local S splitting measurements.

A third swarm, located to the north of the island throughout November 2011, which is likely the result of a collapse of deeper segments of the magma plumbing system (González et al., 2013), provided a large number of larger and deeper (∼20 km) events usable for robust measurements at station CTIG (EH11.EI; Figure 4A). The FPD pattern shows an almost uniform NE-SW pattern, which represents a radial direction outwards from the location of the deeper magma reservoir. This would coincide with the direction of microcracks, which are assumed to be the predominant cause of crustal shear-wave splitting in volcanic regions and which would be oriented radially to overpressurised magma reservoirs (Johnson, 2015; Serrano et al., 2023).

A similar uniformly radial pattern of FPD can be observed during a swarm of high seismicity over 4 days in September 2012 (EH12.I; Figure 4B). Since the events are located at similar depths but further to the south, a similar cause for the swarm can be assumed. Overpressurisation of the deep reservoir is likely to be the reason for the continuing radial alignment of cracks. This is supported by the observation of strong vertical deformation and has been modelled as an intrusive event through joint inversion of GPS and InSAR data (Benito-Saz et al., 2017).

The complete change in pattern a year later (EH13.I; Figure 4C) to a uniformly tangential FPD can be the result of an underpressurisation of certain parts of the reservoir, leading to a change of maximum horizontal compressive stress (Johnson, 2015). This could be caused by a different distribution of island uplift during this intrusion event in comparison to the previous swarm (further to the south and west), even though the crack deformation sources are located in similar places towards the southern tip of the island (Benito-Saz et al., 2017). The slightly shallower hypocentre location of events are likely caused by collapses of shallower segments of the magma plumbing system due to a continuing emptying of the system towards the shallower reservoir. Ongoing subsurface dynamics are evident from the continuing uplift that is observed by GPS measurements throughout 2013. At the same time, events caused by a change of the southern, shallower reservoir around that time were again resulting in shallow incidence angles towards the active stations at the time, thus preventing robust local S splitting measurements.

Afterwards, the general decrease in seismic activity and shear-wave splitting measurements, as well as the lack of strong FPD alignment hint towards the decreased activity in the magmatic plumbing system. A further short intrusive event around March 2014 resulted in nearly 300 additional earthquakes. However, the maximum magnitudes stayed below the limit of 2.5. The relative quiescence is also observed through the stabilisation of vertical GPS measurements and lasts for several years until 2018.

The increase in seismicity, notably of events with a magnitude above 2.5, and shear-wave splitting measurements in 2018 can be caused by a recurrence of activity in the system. Events are predominantly located on- and offshore the southern coast of the island. Measurements are caused by much deeper events (EH18.C; Figures 4D, 6D), hinting towards activity beneath and potential around the deeper magma reservoir. The overall increase of time lag with depth can be caused by multiple anisotropic layers or vertically extended anisotropy due to a broader region experiencing horizontal stress, as well as the presence of melt. This is particularly interesting, since GPS measurements do not show a renewed increase.

The pattern of renewed activity continues in 2019 with events on a broad depth interval (EH19a.C; Figure 4E), but overall showing a less focused FPD. This changes during two seismic swarms in September 2019 (EH19b.C; Figure 4F) and January 2020 (EH20.C; Figure 4G). Events for both swarms are located very deep (>30 km), hinting towards activity in the deeper parts of the magmatic plumbing system. In addition, both swarms are located over a broad area to the west (EH19b.C) and the south (EH20.C) of the island. The FPD patterns are uniformly aligned but show roughly tangential orientation for EH19b.C and roughly radial orientation for EH20.C. These changes can be caused by ongoing dynamic changes in different parts of the plumbing system, consisting of periods of under- and overpressurisation, which do not result in a big overall volume change that would affect the vertical movement of the island.

6.2 La Palma

Since the main seismic swarms prior to and during the eruption were located beneath the central southern part of the island, the location of station EHIG provided a comparably large number of good station-event pairs. In addition, the network was already established before the increase in seismicity in September of 2021 a few days before the eruption, so that shear-wave splitting measurements could be taken from those early events (LP21a.C; Figure 5A). The early shallow swarm (∼10 km) could be used to trace magma rising towards the surface (Carracedo et al., 2022; Cabrera-Pérez et al., 2023). The strong uniform FPD alignment in E-W direction is in general agreement with the predominant FPD from v_P azimuthal anisotropy model of a seismic tomography inversion (Serrano et al., 2023). During the early stages the alignment is likely caused by compressive stress induced by the expansion of the conduit due to the magma influx. During this period the time lag stayed uniformly low at around 0.1 s, which is likely caused by the short event-station distance but can in part be the result of relatively moderate pressures at that time. Tomographic inversion of local seismicity shows an elongated anomaly of high v_P/v_S (>2.0), located in the northwestern part of the Cumbre Vieja (Serrano et al., 2023), which were interpreted as magma reservoirs containing high degrees of partial melt (Kasatkina et al. (2014)), supporting our findings.

Throughout most of the eruption period, the FPD at station EHIG continued to be strongly aligned in E-W direction, though the median event latitude was located slightly further south (Figures 5B–G, 6H). This indicates an ongoing overpressuring of the conduit through continuing magma influx towards the surface. The increase in median time lags (LP21b.ED: 0.17 s; Figure 6L; Table 1) throughout this period also hints towards a stress increase through conduit expansion from the ongoing rising magma. The increase in event depth and median time lag coincides with a change in magma composition due to an activation of a deeper magma source (Plank et al., 2023). During a short period of 5 days, just before a significant increase of seismicity for 1 day on 17/11/2021, the FPD changed to a strong NW-SE alignment. Although this pattern is still likely a result of compressive stress, the deviation from the general pattern indicates an intermediate change in influx, leading to the one-day sharp increase in seismicity, in which the FPD changes back to the predominant E-W alignment.

At the same time, station TBT shows more varying FPD patterns. Because of the small number of results obtained by this station due to its greater distance to the swarms, it is more difficult to derive conclusions from the apparent patterns. However, until mid-November, the FPDs align uniformly in individual time intervals (Figures 5B–D). Since only events from the deep swarm provided good event-station pairs (30–40 km depth), the results focus on deeper structures of the magma plumbing system. The more radial FPD alignment (NW-SE) during the early stages of the eruption hint towards compressive stresses in the deeper structures and the relatively small time lags indicate that most of the longer ray paths to the station do not traverse particularly anisotropic regions (LP21b.ED; Figures 5B, 6K, L). Time lags increase significantly afterwards and FPD alignment changes to E-W (LP21c.ED; Figures 5C, 6K, L) and then tangential NE-SW alignment (LP21d.ED; Figure 5D, 6K). The tangential orientation could be explained by an underpressuring of the deeper conduit due to an emptying of the reservoir during the eruption (e.g., del Fresno et al., 2023), related to a decrease in magmatic influx and subsidence of the island.

After the eruption period, the amount of larger seismicity (M > 2.5) decreases sharply. Lowering the lower magnitude limit to 2.0 provides a higher number of good event-station pairs, but still the overall number of results is too sparse to perform a detailed temporal analysis. Still, a short period of 4 days at the end of March 2022 with an increase of results from the shallow swarm can be observed. There the FPD alignment at station EHIG suggests an increase in compressive stress around a depth of 10 km. However, the total number of results is low and generally consists of small magnitude events (M < 2.5). In addition, at station TBT only one result can be observed, which comes from an event shallower than any shear-wave splitting observations during the eruption (<30 km) and in general no results from deeper events can be observed after the end of the eruption. This in combination with the overall lack of deeper seismicity is a strong hint that there are no significant stress changes around the deeper part of the magmatic plumbing system, likely caused by a stop of deeper magmatic movement.

6.3 Anisotropy and volcanic unrest

The time lags observed during and after the eruptions on both islands are larger than normalised values for most crustal rocks, indicative of the higher heat flow in volcanic environments (Crampin, 1999; Bianco et al., 2006). Temporal variations can be caused by changes in pore pressure or in crack orientations (Miller and Savage, 2001). We observe high scatter of time lag values, resulting in average variations that are less significant between individual time intervals than the changes in FPD. This phenomenon is caused by the δ t dependency on both the path length as well as the amount of anisotropy on the path, making it a less robust indicator of subsurface stress regime variations (Miller and Savage, 2001, and references therein).

In both islands, we can observe changes in FPD. They are more pronounced in El Hierro, which can be explained by the longer observation time, spanning various intrusive events and island inflation, as well as intervals of relative quiescence, but are observable even in La Palma over the shorter observation time. This includes the ∼60°-flip between EH12.I and EH13.I on El Hierro (Figures 4B, C), as well as the ∼40°-flip between LP21c.ED and LP21d.ED on La Palma (Figures 5C, D). Similar complete flips in FPD of up to 80° have been observed at Mount Ruapehu in New Zealand, being attributed to repeated filling and depressurisation of the magmatic dyke system leading to changes in the local stress field superimposing on the regional stress field (Gerst and Savage, 2004). Likewise, flips of up to 90° were observed around eruptions of Mount Etna in Italy and attributed to realigned microcrack distributions due to an increase of pore-fluid pressures to critical levels (Bianco et al., 2006). FPD flips of ∼70°–90° could also be observed in Montserrat, showing local stress field changes due to conduit pressurisation (Roman et al., 2011).

As discussed, our preferred interpretation for changes in splitting parameters (especially FPD) is due to temporal variations, but especially in El Hierro spatial variations cannot be ruled out due to significant changes in event coordinates between different swarms. However, even though the swarm origins for EH12.I and EH13.I are at slightly different depths and have marginally different ray paths, an interpretation of spatial variation seems unreasonable over such small length-scales. A similar conclusion was drawn in a study around the 2001 eruption of Etna, Italy (Bianco et al., 2006).

Previous studies often observe significant pre-eruption changes in splitting parameters (Miller and Savage, 2001; Gerst and Savage, 2004; Bianco et al., 2006) emphasising the need of active dense seismic networks around active volcanoes. While our observations focus more on the time during and after the eruption, the discrepancy between available seismicity and local splitting results on El Hierro during the initial stages of the eruption (Figure 7A) support that idea. This is of particular importance, since it contributes to a better understanding of the changes in different parts of the plumbing system, which do not necessarily result in big overall volume changes that would affect inflation or deflation around a volcanic region, or can be registered by other visible features on the surface, such as outgassing or erupted lava.

7 Conclusion

This study shows the dynamic changes of the magmatic plumbing system before, during and after a volcanic eruption through the measurement of seismic anisotropy patterns. We observe clear indications of intrusive events in the 2 years following the eruption on El Hierro, as well as the week preceding the eruption on La Palma. After an interval of relative quiescence, an increase of results from El Hierro starting in 2018 hint towards renewed subsurface activity within the deeper parts of the system, which affected the rate of overall seismicity but did not show significant changes in vertical GPS measurements at the same time. Activity on La Palma following the eruption has been low over a period of one and a half years, but an activity increase akin to patterns observed in El Hierro is not unlikely.

The measurement of seismic anisotropy provides useful additional constraints that can be combined with observations from other geophysical and geochemical methods. Limitations arise due to requirements of minimum event magnitude and steep event-station ray paths. Ideally, seismic networks should be operational before the volcanic eruption to be able to investigate preceding seismic swarms indicative of activity within the system. Although a detailed temporal analysis is not as straightforward to interpret as, for example, variations in local seismicity or the GPS measurements of vertical movement, it serves as an important addition as it sheds light on subsurface changes not visible by other methods. Ultimately, it helps improve our understanding of the stress field around the plumbing system and adds robustness to models that attempt to forecast the likelihood of an eruption and/or prolonged significant seismicity following seismic crises in regions subjected to intrusive and eruptive activity.

Data availability statement

Publicly available datasets were analyzed in this study. They were provided by Instituto Geografico Nacional (DOI: 10.7914/SN/ES). This data can be found here: https://www.ign.es/web/ign/portal/sis-area-sismicidad.

Author contributions

DS: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. GS: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Writing–review and editing, Data curation. RR: Resources, Validation, Writing–review and editing. PG: Resources, Validation, Writing–review and editing. RA: Resources, Validation, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work is a contribution to projects SIGHT (Ref. PTDC/CTA-GEF/30264/2017), GEMMA (Ref. PTDC/CTA-GEO/2083/2021), and RESTLESS (DOI: 10.54499/PTDC/CTA-GEF/6674/2020). This work was also funded by the Portuguese Fundação para a Ciência e a Tecnologia (FCT) I.P./MCTES through national funds (PIDDAC)—UIDB/50019/2020 (https://doi.org/10.54499/UIDB/50019/2020), UIDP/50019/2020 (https://doi.org/10.54499/UIDP/50019/2020) and LA/P/0068/2020 (https://doi.org/10.54499/LA/P/0068/2020). We thank Spanish Agencia Estatal de Investigación projects COMPACT (Ref. PID2019-104571RA-I00) and Volta-Motion (Ref. PID2022-139159NB-I00) funded by MCIN/AEI/10.13039/501100011033 and “FEDER Una manera de hacer Europa”. Research activities of the CSIC staff during the eruption were funded by CSIC -CSIC-PIE project PIE20223PAL008.

Acknowledgments

We would like to thank the editors and the two reviewers for their useful suggestions to improve the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1299338/full#supplementary-material

SUPPLEMENTARY TABLE S1 | Details of seismic stations used in this study.

SUPPLEMENTARY TABLE S2 | Details of shear wave splitting results for El Hierro and La Palma. Note that the results for El Hierro were first presented in Schlaphorst et al. (2023), but will be listed here again for the completeness of this paper.

References